Overlapping speckle correlation algorithm for high-resolution imaging and tracking of objects in unknown scattering media

Optical imaging in scattering media is important to many fields but remains challenging. Recent methods have focused on imaging through thin scattering layers or thicker scattering media with prior knowledge of the sample, but this still limits practical applications. Here, we report an imaging method named ‘speckle kinetography’ that enables high-resolution imaging in unknown scattering media with thicknesses up to about 6 transport mean free paths. Speckle kinetography non-invasively records a series of incoherent speckle images accompanied by object motion and the inherently retained object information is extracted through an overlapping speckle correlation algorithm to construct the object’s autocorrelation for imaging. Under single-colour light-emitting diode, white light, and fluorescence illumination, we experimentally demonstrate 1 μm resolution imaging and tracking of objects moving in scattering samples, while reducing the requirements for prior knowledge. We anticipate this method will enable imaging in currently inaccessible scenarios.

Supplementary Fig. 2 Algorithm flow chart of the overlapping speckle correlations.Two of the detected N speckle images Ii and Ij are numerically separated into envelopes Ei and Ej and speckles Si and Sj.The Hadamard product is performed on Si and Sj to obtain the overlapping speckle, and all pixel values of the overlapping speckle are summed to obtain a value Cij of the object's autocorrelation.The cross-correlation  is performed on Ei and Ej to obtain the position (xij, yij) of value Cij on the object's autocorrelation.Combining the calculated value and its position information, a pixel Cij(xij, yij) of the object's autocorrelation is constructed.The object's autocorrelation Cij(xij, yij) is completely constructed by successively changing i and j from 1 to N. The object image O is recovered from the autocorrelation through a phase-retrieval algorithm.

Supplementary Note 5: Assumption of uniform IS(x, y)
To confirm the feasibility of the assumption above Eq.( 5) that the intensity distribution IS(x, y) is nearly uniform within a small region containing the object, we detected the intensity distributions on the object plane by removing the latter half of the parafilm.
As shown in Supplementary Figs.

Supplementary Note 6: The imaging depth limit
To explore the imaging depth limit of speckle kinetography, imaging of an object consisting of three 25 m width lines moving inside parafilm of 2 to 30 layers is experimentally implemented.The structural similarity index measurement (SSIM) of the constructed autocorrelation decreases as the sample thickness L increases (Supplementary Fig. 6a).When the sample thickness reaches 28 layers, the object's autocorrelation becomes too blurry to image.Therefore, the maximum imaging depth here is about 26 layers of parafilm, corresponding to 19.9ls and 6.2lt.To determine the main limiting factors, the envelope expansion, the relative displacement error and the speckle contrast versus the sample thickness are analyzed.The full width at halfmaximum W of the envelope broadens as the sample thickness L increases (Supplementary Fig. 6b).But the entire envelope can still be measured when the parafilm is 28 layers, so it does not have obvious impact on imaging.The relative displacement error E fluctuates slightly as L increases (Supplementary Fig. 6c).But the fluctuation keeps within ±4 camera pixels, which is smaller than the pixel size of the object's autocorrelation, corresponding to 16 camera pixels in these experiments.Thus, it does not cause misplacement to autocorrelation pixels.The speckle contrast decreases as L increases and it quickly drops below 0.1 from 26 to 30 layers (Supplementary Fig. 6d).The low contrast causes the speckles to overlap with each other even when the corresponding objects do not overlap.In this case, the small values of the object's autocorrelation become large, which is consistent with the blurred autocorrelation constructed at 28 layers (Supplementary Fig. 6a).It causes the imaging to fail.In conclusion, although the above three factors all limit the imaging depth, the speckle contrast dominates.To eliminate the anisotropic effects of the scattering sample, imaging of the abovementioned 25 m widths object embedded in homogeneous static isotropic polyethylene foams (Supplementary Figs.7a and 7b) is experimentally performed.
When the polyethylene foam is 15 mm thick, the speckle images with entire envelope and sufficient speckle contrast are recorded under narrowband illumination (singlecolour LED of 625 nm nominal wavelength and 17 nm bandwidth, 920 mW) and broadband illumination (research arc lamp source of 260 to 2500 nm wavelength, 150 W), as shown in Supplementary Figs.7c and 7d.The constructed autocorrelations and recovered object images are shown in Supplementary Figs.7e-7h.However, when the polyethylene foam is 20 mm thick, the contrast of the speckle images is too low to image.Therefore, the maximum imaging depth here is between 15 mm to 20 mm where the C(x, y) in the 3 rd step is the measured object's autocorrelation.
The input for the next (i+1) iteration, oi+1(x, y), is obtained from the output of the i th iteration, oi ' (x, y), by imposing physical constraints on the object image, which is real and non-negative in our implementation.Two types of implementations of these constraints are used in the algorithm [9][10][11] , termed the Hybrid Input-Output (HIO) algorithm: where Γ is the set of all points (x, y) on oi ' (x, y) that violate the physical constraints, and β is a feedback parameter that control the convergence properties of the algorithm.Following Katz et al. 9 , a few thousand iterations of the HIO algorithm were ran with a decreasing beta factor from β = 2 to β = 0, in steps of 0.04.For each β value, 40 iterations of the algorithm were performed.The result of the HIO algorithm was fed as an input to additional 40 iterations of the Error reduction algorithm to obtain final result.
To assure faithful reconstruction of each image, several different runs of the algorithm (from 10 up to 60, typically 20) were performed with different random initial conditions, and the reconstruction having the closest Fourier spectrum to the measured autocorrelation Fourier transform (lowest mean-square-error) was chosen as the final reconstructed result.

Supplementary Note 3 : 4
Autocorrelation construction with a simplified trajectory Two examples of autocorrelation construction with trajectories of U shape and T shape are shown in Supplementary Fig. 3.The construction method is described in the Methods of the manuscript.Supplementary Fig. 3 Autocorrelation construction with a simple trajectory.The object image can be reconstructed as long as the autocorrelation of the trajectory is large enough to cover the object's autocorrelation with sufficient resolution.a Speckle image from a hidden object.The trajectory extracted from the cross-correlations between every two envelopes is in a U shape.The overall size of the trajectory is slightly larger than the speckle pattern.b Autocorrelation of the U-shaped trajectory in (a).c The constructed autocorrelation of the object with a sampling area determined by (b).d The object image recovered from (c).e Speckle image from a spontaneously moving object.The trajectory is recovered from the crosscorrelations between every two envelopes.f Autocorrelation of the partially selected trajectories shown as a blue dotted T shape in (e).g The object's autocorrelation with a sampling area determined by (f).h The object image recovered from (g).Scale bars: 50 m.Recovered trajectories in the experiments.The relative positions are the white spots in the centre of the yellow squares.The red solid lines represent the real trajectories.The width of the red line is set approximately equal to the pixel width of the constructed autocorrelations.a The recovered trajectory of the object consisting of two 10 m width lines and moved in 22 layers of parafilm.b The recovered trajectory of the object consisting of three 5 m width lines and moved in 14 layers of parafilm.c The recovered trajectory of the object consisting of six 1 m width lines and moved in 6 layers of parafilm.d The recovered trajectory of a number-shaped 5 and moved in 14 layers of parafilm.e The recovered trajectory of the fluorescent beads.Scale bars: 10 m.
5a-5c, a transmissive object consisting of three lines was illuminated by a narrowband halo with intensity distribution IS(x, y) emitted from parafilm of 3, 7 and 11 layers, respectively.Supplementary Figures5d-5fshow the number-shaped object illuminated by a broadband halo with intensity distribution IS(x, y) emitted from parafilm of 3, 7 and 11 layers.Although the intensity distributions shown in Supplementary Fig.5are not completely uniform, the intensity fluctuation is minimal relative to the overall intensity.Thus, the intensity distribution IS(x, y) within a small region containing the object can be assumed to be nearly uniform.Supplementary Fig. 5 Intensity distribution on the object plane.a-c Intensity distributions of an object consisting of three 5 m width lines that are illuminated by a narrowband halo with intensity distribution IS(x, y) emitted from parafilm of 3 (a), 7 (b) and 11 (c) layers.d-f Intensity distributions of the number-shaped object illuminated via a broadband halo with intensity distribution IS(x, y) emitted from parafilm of 3 (d), 7 (e) and 11 (f) layers.Scale bars: 10 m.

Supplementary Fig. 6
The influence of medium thickness on four factors.Imaging of an object consisting of three 25 m width lines moved in parafilm of 2 to 30 layers is experimentally performed.a The SSIM of the constructed autocorrelation versus the sample thickness L. The insets are the constructed autocorrelations from 2 to 26 layers with an interval of 8 layers.Their corresponding SSIMs are respectively marked with red dots.The autocorrelation (the last inset) becomes too blurry to image when L reaches 28 layers.The corresponding SSIM is obviously smaller than other SSIMs and is marked with a red dot.b The FWHM W of the envelopes versus L. c The relative displacement accuracy E versus L. d The speckle contrast versus L. The speckle contrast is below 0.1 when L reaches 28 layers.
foams, corresponding to 5.1~6.8ls and 5.0~6.7lt.Notably, the maximum imaging depths in polyethylene foam (5.0~6.7lt) and parafilm (6.2lt) are consistent after eliminating the anisotropic effects.Therefore, we conclude that the imaging depth limit of speckle kinetography is about 6lt.The scattering properties of the scattering samples are described in Supplementary Note 10.Supplementary Fig. 7 Imaging of an object embedded in 15 mm thick polyethylene foams.a and b Front (a) and side (b) photographs of the polyethylene foams with 15 mm thickness.c One of the recorded speckle images under narrowband illumination.d One of the recorded speckle images under broadband illumination.e The object's autocorrelation constructed from the recorded speckle images in (c).f The object's autocorrelation constructed from the recorded speckle images in (d).g The object image recovered from (e).h The object image recovered from (f).Scale bars: 50 m.